arXiv:1501.01662vl [physics.plasm-ph] 7 Jan 2015 


Approximate semi-analytical solutions for the 
steady-state expansion of a contactor plasma 


E. Camporeale 1 , E. A. Hogan 2 , E. A. MacDonald 3 

1 Center for Mathematics and Computer Science (CWI), Amsterdam, Netherlands. 

2 University of Colorado, Boulder, CO, USA 

3 NASA Goddard Space Flight Center, Greenbelt, MD, United States 
E-mail: e.caLmporeale@cwi.nl 

Abstract. We study the steady-state expansion of a collisionless, electrostatic, quasi¬ 
neutral plasma plume into vacuum, with a fluid model. We analyze approximate 
semi-analytical solutions, that can be used in lieu of much more expensive numerical 
solutions. In particular, we focus on the earlier studies presented in Parks and Katz 
(1979) Jlj, Korsun and Tverdokhlebova (1997) [2], and Ashkenazy and Fruchtman 
(2001) |3]. By calculating the error with respect to the numerical solution, we can 
judge the range of validity for each solution. Moreover, we introduce a generalization 
of earlier models that has a wider range of applicability, in terms of plasma injection 
profiles. We conclude by showing a straightforward way to extend the discussed 
solutions to the case of a plasma plume injected with non-null azimuthal velocity. 
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1. Introduction 

The modeling of the expansion of a plasma plnme in the vicinity of a spacecraft has been 
intensively investigated in the last several years. In particular, in the context of space 
electric propulsion systems, such as ion or Hall thrusters, a correct characterization 
of the emitted plasma plume is of crucial importance to avoid interactions between 
the energetic particles and spacecraft surface, and thus to potentially prevent severe 
damages to the spacecraft pQElEltZ!. Similarly, in plasma contactor technology (for 
instance, electrodynamics tether applications), it is important to predict the shape and 
geometry of the plasma plume emitted by a spacecraft mm- 

Due to the large difference in density between the plasma plume and the background 
magnetospheric plasma, the physics of the plume expansion is often well described by 
the expansion of plasma in vacuum, which has also been investigated thoroughly for 
several different geometries [10 4 HU- The first-principle numerical calculation of the 
steady-state profiles reached by the expansion of a plasma plume in a large domain 
is very challenging, due to the wide separation of scales involved. For instance, the 
plasma Debye length can vary by several orders of magnitudes in few hundred meters 
along the expansion trajectory. Therefore, numerical simulations are usually constrained 
to either the near- or the far-held expansion, i.e., to a restricted simulation domain 

A very appealing alternative to expensive numerical simulations is, of course, to seek 
for analytical solutions of the plume expansion. This is often done by assuming that 
the plasma obeys simplified fluid equations E El El ESI EDI EH- 

In this paper we analyze and discuss the self-similar solutions presented in Refs. EM. 

It is important to notice that the self-similar approach used in such models is different 
from the standard procedure followed in similarity methods described, for instance, in 
classical textbooks m In fact, the standard approach is based on the idea of reducing 
the dimensionality of a set of partial differential equations, through an appropriate 
variable transformation, that can be identified by exploiting the symmetries of the 
system. On the other hand, the models described in this paper retain their full- 
dimensionality, and are based on a variable transformation constructed on the principle 
that one variable is constant along fluid streamlines. Another crucial point to notice, 
is that the solutions presented in Refs. E13 E] are approximate solutions. As we will 
show, the self-similarity assumption is indeed inconsistent with the fluid model employed 
and, as such, the solutions obtained do not exactly solve the initial set of equations. 
However, such solutions, albeit inexact, have the advantage that they can be quickly 
evaluated on an arbitrary large domain. For this reason, they can still be valuable, if 
one is able to estimate how large their error is with respect to the true solution, and, of 
course, if such error is reasonably small. 

The aim of this paper is twofold. First, we analyze the earlier solutions of Refs. EM 
and we measure their errors with respect to the numerical solution of the underlying 
fluid equations. Second, we present a generalization of such solutions, in an attempt to 
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provide a more flexible family of solutions that can have a wider range of applicability. 
The paper is organized as follows. Section 1 presents the mathematical model, the 
self-similarity assumption, and the resulting set of equations. In Section 2, we describe 
the approximate solutions presented in Refs. mim and we comment on their range 
of validity and respective errors. In Section 3, we describe a new class of solutions, 
that is based on a generalization of previous models, and we show that they are indeed 
applicable to a wider class of situations. In Section 4 we show how to easily extend the 
solutions that were derived for a plasma with no azimuthal velocity to a more general 
case. Finally, we draw conclusions of this study in Section 5. 

2. Mathematical model 

We study a system composed by collisionless, singly charged, ions and electrons, in 
a steady-state, axisymmetric configuration. We assume that electrons inertia can be 
neglected, quasi-neutrality holds, the plasma is electrostatic, and there is no background 
magnetic field. For simplicity, we assume cold ions although a polytropic equation of 
state for ions might be easily incorporated in the model. By employing cylindrical 
coordinates (r, z, 6) the ion continuity equation and the conservation of ion and electron 
momentum read: 
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where the symmetry condition ^ = 0 has been used, n is the density, u z and u r are 
the ion axial and radial velocities, <f is the electrostatic potential, and p e is the electron 
pressure. Eqs. m-m are complemented by a polytropic equation of state for the 
electrons: 

Pe = n 7 (7) 

with 7 the polytropic index. Quantities have been normalized as follows: velocities 
to a/tt7 mi (with To a reference electron temperature and m* the ion mass), the 
electrostatic potential to T 0 /e (with e the elementary charge), density to a reference 
density n 0 , lengths to a characteristic length R, and pressure to n 0 T 0 . In general, the 
value of the polytropic index 7 depends on the degree of ionization of the plasma [23]. 
For instance, [7] have reported experimental results for thrusters, with 7 lower than 
5/3, and [24] have used the value 1.3 for their PIC simulations. However, as shown in 
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Ref. |23j. 7 tends to the neutral gas theoretical value for a fully ionized plasma, which, 
for simplicity, is the case treated in this paper. Hence, for all the cases presented here 
7 = 5/3, i.e. the adiabatic constant of the monoatomic perfect gas. 

The model in Eqs. 0-10 must be interpreted as a boundary value problem in the two- 
dimensional plane (r,z). Since the equations involved are first order partial differential 
equations, one is allowed to specify the boundary conditions at r = 0 and z — 0 , which 
determine the solution in the whole domain. Also, note that the electron pressure and 
the electric potential can be substituted in Eqs. (J2])-(EJ) , and post-processed after the 
solution for n, u r , u zi and ug has been obtained. 


2.1. Self-similar solution 

We now elucidate the procedure to seek for a self-similar solution of Eqs.(Q])-(j5]), 
extending the derivation in [25] to the case Ug ^ 0. First, we introduce a change of 
variable (r, z) —> ( r],z ), with p = -A A, and a(z ) an unspecified function. Then, we 
assume that the unknowns n,u z ,u r ,ug can be factorized as: 

n(p,z) =n c (z)n t (r /) ( 8 ) 

u z (jji z ) — u c (z)u t (p) (9) 

u r (v, z) = pa'(z)u z (p, z), ( 10 ) 

( 11 ) 


T) 

M r ,v) = 

a(z) 

where prime indicates differentiation, and E is an arbitrary constant. The definition 
for u r , Equation (TTOT) . follows from assuming that lines of constant p correspond to 
streamlines, that is 
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As we anticipated, the set of equations ([I])-© is actually inconsistent with the separation 
of variables assumed in (l 8 l)- (fTTTh That is, the only solution that satisfies such separation 
of variables is the trivial solution with n, u z ,u r , Ug all constant. However, the separation 
of variables is still worth considering if one can derive a class of approximate solutions 
for our model that are much faster to compute than the exact numerical solution, and 
its error relative to the full solution is small. 

It is straightforward to prove that Eq.Q is automatically satisfied by using the 
factorization assumption in Eqs. ©-CD- Substituting Eqs.(l8])- (fT0|) into Eqs.(pQ)-(j5]) and 
separating each equation in p and z dependent terms, one can get the following set of 
equations for u c (z), n c (z), and a(z): 
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and the the following equations for n t {rj) and u t (rj): 

2-D(7-l) _ —2D 
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where A, B , C and D are arbitrary (separation) constants. Details of the calculation 
can be found in the Appendix. Note that Eq. (]T3]) descends from the continuity equation 
CO), Eqs. (HU), (JED, (ITT]) from the momentum equation in the z direction (J3]), and Eqs. 
(1T6]) and (jT8T) from the momentum equation in the r direction (]2]h Clearly, we have a 
set of 6 equations for only 5 unknowns, and the system is overdetermined. 


3. Approximate solutions for ug = 0 


In order to derive an approximate solution of the model, a common approach is to 
assume u c = const and to prescribe a certain profile for ii t . In this way one can still 
satisfy some (but not all) of the equations (TT3|) - (jT8|) . For simplicity, and to be adherent 
to earlier literature, we consider in this section the case where E — 0, and thus ug = 0. 
All of the solutions discussed here satisfy the continuity equation and the conservation 
of momentum in r, but not in 2 . From equations (113]) . (113]) . and (ITS]) . it follows that 
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In addition to the value of u c and the profile of u t , Eqs. ffTTJ])- (12T1) require to specify 
the value of C, a(0), a'(0), and n c (0). As we will see some of these free parameters 
depend on the boundary conditions at z = 0, i.e., the boundary from which the plume 
is injected. We assume the following boundary conditions: 
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Eqs. (j 22]) and (124]) are simply normalization constraints. Eqs. ([23]) requires that the 
prohle of the density at injection decreases by a factor of 100 within the characteristic 
length R. Eqs. (1251) - (127]) are symmetry conditions at r = 0. 

In this section, we describe the solutions presented in m. 0.0. and we propose a new 
class of solutions that include earlier models as limiting cases and produce an higher 
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degree of flexibility in terms of plume injection profiles that can be represented at z = 0. 
Such approximate solutions provide the density and velocity profiles at the boundary 
z — 0 which are used to calculate the (exact) numerical solution. The latter is evaluated 
on a domain [r, z] = [0, R] x [0,80], with R = 50 and grid size Ar = A z = 0.2. The 
derivatives in r and z are discretized with a fourth order central (5-points stencil) and 
a third order upwind (4-points stencil) difference scheme, respectively. 

Since we are only interested in the steady state solution, we have used the following 
marching scheme for solving numerically Eqs. (fT|) -(l5j). First, the solution is calculated 
on the partial domain [0, R] x [0,4Ax:] with an iterative Newton-GMRES solver [26]. 
The solution is then extended to the full domain by adding one row of cells in z at the 
time and calculating the new solution at each step. This is possible since we use an 
upwind discretization in z that requires only the points where the solution has already 
been calculated. 

For each approximate solution we evaluate the error with respect to the full numerical 
solution. We define the following two measures of error: 
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which are the L\ norm of the relative error of the radial and axial fluxes, where the 
superscripts a and n indicate the approximate and the numerical solutions, respectively. 
Each of the solutions presented in Refs. dm G] makes a different assumption on ut, 
and we proceed to discuss them separately. 


3.1. Parks and Katz (PK) solution 


The solution proposed in [1] assumes a constant u t 
profile: 
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where we have used n c (0) = 1, imposed by Eq. (l22lh The boundary condition (f24l) 
allows to solve for G and Eq. (1241) yields the relationship 

a(0) = u c a'( 0) 


(34) 
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In summary, the only two free parameters are u c and a'(0) and the profile for a(z) can 
be calculated by solving numerically Eq. fl20|) . Note that the profiles of u r and n at the 
injection boundary z = 0 do not depend on the particular choice of u c and a'(0). Figure 
ID shows the profiles for the density and the radial velocity at z = 0 as a function of r. 
Clearly, the radial velocity increases linearly with r. An important feature of the PK 
solution is that, for any given value of z, the solution is constrained to a certain range 
in r. Indeed, one can see from Eq. (l3T]) that the density becomes a complex quantity 
(i.e. unphysical) for a large enough r. In Figure [H we have artificially set n = 0 in 
the regions where the density becomes imaginary. Figure [2] shows the density profile in 
the (z, r) plane, for different values of the axial velocity u z = u c , in logarithmic scale, 
for a fixed value of a'(0) = 0.2. Once again, notice that the density profile at z = 0 
is not a function of u c . The errors e r and e z as a function of u c are shown in Figure 
|3] For u c > 20, e r and are smaller than 1%. Note, however, that here and in the 
next similar figures, the error has been calculated in the restricted region |r| < R. The 
region where the density is artificially set equal to zero is excluded by the evaluation of 
the error. The reason is that the marching routine described above is not very robust 
when a sharp gradient in the density is included as boundary condition and, therefore, 
it is not straightforward to obtain a numerical solution that extends to |r| > R. 


3.2. Ashkenazy and Fruchtman (AF) solution 


The AF solution J3] assumes a conical velocity profile at the injection, such that 

a, — (1 - a'(O)V) 1 


(35) 


The profile for n t is given by solving Eq. (j2Tj) : 

nt = ( c ° - c Zm log ^ + (a ' (0))V) ) 7 ’ (36) 

where Co is an integration constant that we set equal to 1. The boundary condition 
(l22j) imposes n c ( 0) = 1, while Eq. fl23jl yields an equation that must be solved to obtain 
the value of C. The remaining free parameters u c , a(0), a'(0) obey, via the boundary 
condition (J24J), the relationship: 
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For simplicity, we fix a'(0) = 0.2, allowing u c to vary, and calculate a(0) through Eq. 
m- As usual, a(z) is calculated via Eq. fl20)) . In summary, the prohles for n, u r , and 
u z are: 


n (r, z) 
u r (r,z) 


o(0) 

a(z) 


1 -AcAiog |1+ 


a(z) 


a\z)u c 


1 - 


2 ( a '( 0)) 2 

g(7 ~ 1) 

2 ( a '( 0)) 2 


log 1 + 


a'( 0 ) , 

a(z) ‘ 

a'(0) 
a(z ) 


7-1 


7-1 


(38) 

(39) 
















Approximate semi-analytical solutions for the steady-state expansion of a contactor plasma^ 


u e (r, z ) 



(40) 


Figure 0] shows the density profile for the AF solution at the injection boundary z — 0, 
as function of r, for the values u c = 2, 10, 20, 100. Similarly to the previous model, the 
density becomes imaginary for r larger than certain values. The corresponding profiles 
for u z and u r are plotted in Figures [5] and [61 respectively. Figure [7] shows the two- 
dimensional density profile in the (z, r) plane, in logarithmic scale for four values of u c . 
One can notice a qualitative similarity with the profiles of the PK solution, especially 
for large u c . The errors for the AF solution are shown in Figure [Eland are comparable 
with those obtained for the PK solution (Figure E]), he., the errors are smaller than 1%, 
for u c > 20. 


3.3. Korsun and Tverdokhlebova (KT) solution 


The solution proposed in [2] assumes the profile 
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Moreover, by setting D = —2, the KT solution satisfies Equation (ITT)) . By applying the 
boundary condition (I T2\i and Eq. (fT51) . it follows that a(0) = n c (0) = 1. In summary, the 
KT solution not only satisfies the equations of continuity and momentum in r, but also 
Eqs. (jT5|) and (TTTj) . i.e., it satisfies 2 out of 3 equations for the conservation of momentum 
in 0 direction. The profiles for n, u r , and u z are: 
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The boundary condition (fTTl) provides the expression for C = 2 (l — / R 2 = 

—0.0792, and a(z) is calculated via Eq. (l20]h Hence, the only free parameter is u c which 
is related to a'(0) through the relationship 
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that is obtained by applying Eq. fl24]h 

Similarly to the PK solution, the density and radial velocity profiles at the injection 
boundary z = 0 are not a function of u c . They are shown in Figure [9] The profile for 
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u z is a linear function of u c . Three examples for the values u c = 1, 10, 100 are shown in 
Figure flOl Four examples of two-dimensional density profiles are shown in logarithmic 
scale in Figure fill for u c = 2, 10, 20, 100. The errors e r and as a function of u c are 
shown in Figure [12] Notice that they are much higher than the respective errors for the 
PK and AF solutions. This might be surprising since the KT solution satisfies equations 
(Eli and Eli, while the PK and AF solutions do not. However, given the nonlinear 
nature of the model, one cannot predict how large the errors will be with respect to the 
true solution. Moreover, as discussed in the Appendix, Eqs. (TT5l) and (fT71) follow from 
the assumption that Eq. (TT4l) holds. Hence, satisfying equations (jT5|) and (TT7l) does not 
automatically yields a lower error, since Eq. (TT4li is still not satisfied. 


3.f. Discussion on the PK, AF, and KT solutions 


We have noted a certain similarity between the PK and AF solutions for large u c . This 
is not surprising if one realizes that the AF solution tends exactly to the PK solution 
in the limit u c —> oo. Indeed, in this limit, Eq. (1371) requires that a(0) —> oo and/or 
a'(0) —> 0. Taylor expanding the logarithmic term in n and u r [Eqs. (138|) - (139|) ] one has 


(47) 



and Eqs. (I38lb(j40l) reduce exactly to Eqs. (l3TTh(j33lh 

Although the PK and AF solutions yield relatively good results in terms of small errors 
with respect to the numerical solutions (especially for large u c ), they are severely limited 
in the choice of density and velocity profiles at the injection. In fact, the value of u c 
determines, by construction, the injection profiles that, in turn, determine the region 
in which the plume propagate. From Figures [2] and [3 it is evident that the cone of 
propagation becomes more and more collimated as u c increases. In Figure [I3l we show 
the angle (in degrees) of the isocontour n = 0 with respect to the axis r = 0, as 
a function of u c . Note that the cone of propagation is identical for the AF and PK 
solutions. Also, although the isocontour n — 0 is not exactly a straight line in the ( z , r) 
plane, it can be considered approximately straight in the region plotted in Figures [2] 
and 0 We emphasize that this cone of propagation, and the relative angles shown in 
Figure [[3] are an intrinsic characteristic of the PK and AF solutions, which depends 
solely on the value of u c . On the other hand, the user of such approximate solutions 
might need to be able to choose a different injection prohle than the one imposed by 
these solutions, since these depend on the source generating the contactor plasma. In 
particular, as highlighted in Figure [131 large injection angles are not possible for a large 
injection velocity, for the PK and AF solutions. Differently from PK and AF, the KT 
solution does not involve a propagation cone. Indeed, from Equation (l43]h it is easy 
to see that the density tends to zero only asymptotically for r —)■ oo. Unfortunately, 
although the KT profiles offer an alternative to the PK and AF solutions, they yield a 
much larger error and therefore are not preferable. 

In the following section we introduce a new class of solutions that allows a rather 
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wide choice of injection profiles, still yielding relatively small errors with respect to 
the numerical solutions. 


4. A new class of solutions 


The new class of solutions is based on a generalization of the KT and PK solutions. 
Interestingly, their error is much smaller than the errors yield for the KT solution, and 
comparable to the ones of the PK and AF solutions. As in previous models, we still 
assume a constant u c . We also assume the following profiles for n t and up. 
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Note that the KT solution belongs to this class of solution, for the particular case 
D = —2 and F = 1, while for D — 2/(y — 1) and F — 1, one recovers the PK solution. 
Also, one can verify that Eqs. (l48l) -(l49l) satisfy Eq. (fTTlh 
The solution reads: 
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It follows that four parameters can be freely chosen between D , F, a(0), a'(0), and u c . 
In order to simplify the study and to be consistent with the previous models, we fix 
a'(0) = 0.2, F — 1, and we vary D and u c (solving Eqs. (l54li - (l55li for a(0) and C ). Note 
that the profiles of n, u r , and u z at the injection boundary z = 0 do not depend on 
a'(0). 

Figure [H] shows the density profile at £ = 0 for D = —7 (green), -3 (black), 1 (red), 
5 (blue), as a function of r. Clearly, the value of D controls the width of the injected 
density, that can vary from Gaussian-like to a profile closer to a step-function. Also, 
similar to the KT solution, the density tends to zero asymptotically when D < 0, for 
r —> oo. In this case, however, the value of D determines the asymptotic scaling. On 
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the other hand, for D > 0, one has imaginary values that must be artificially set equal 
to zero. Figures [J5] and |T 6 ] show the profiles of u r and u z , respectively, for D = —7 
(green), -3 (black), 1 (red), 5 (blue), obtained with u c = 2. Note that whilst the profiles 
of n and u r do not change by varying u c , the axial velocity u z is a linear function of 
u c , i.e., the profiles in Figure [T 6 ] are linearly rescaled by changing the value of u c . In 
Figure [17] we present a representative example of the two-dimensional density profile, in 
logarithmic scale, for u c = 2, 10, 20, 100, and D = —5. In order to show the versatility 
of this class of solutions, we show more examples of the two-dimensional density profile 
in Figure [TB1 for u c = 20 and varying I) = —7, -3, 1 , 5. Clearly, the cases with D < 0 
do not present a propagation cone. 

Figures [T9l and 1201 present the errors e z and e r , as a function of u c , for D = —7 (blue), 

-3 (green), 5 (red). As anticipated, the errors are comparable to the ones for the PK 
and AF solutions, and are monotonically decreasing with increasing u c . 


5. Extension to the case ug 0 


By still considering a constant u c , Eq. 


is modified as: 
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Clearly, from Eq. flTTl) . the value of E determines the magnitude of the azimuthal velocity. 
We plot in Figure [21] the ratio between the total velocity v = (vf + u 2 + uff) 1 ^ 2 and the 
azimuthal velocity uq, at the injection boundary z = 0 , as a function of r, for different 
values of E. Note that ug increases linearly with r, and v/ug does not depend on D. 
The inclusion of a non-null azimuthal velocity at injection causes a distortion of the 
density profile. Figure [22] shows the two-dimensional density profile in (z,r), for E = 0, 
500, 1000, 2000, D = —10, and F — 1. The notable effect of the azimuthal velocity 
is that the density decreases much more sharply along the axial direction, when ug is 
sufficiently large. This effect is even more evident when D > 0, i.e., when the plume is 
bounded by a propagation cone. Figure [221 shows the two-dimensional density profile in 
(z, r) for D = 3 and F = 1 (i.e. for the PK solution). Clearly, the cone of propagation 
is distorted and the isocontour n — 0 becomes more and more curved, with increasing 
E. 


6. Conclusions 

We have analyzed in detail the self-similar solutions presented in ESS], that represent 
an approximation solution for the fluid modeling of the steady-state, axisymmetric 
expansion of an electrostatic, quasi-neutral, collisionless plasma plume. It is important 
to remind that the semi-analytical solution can be computed in a negligible fraction of 
the time needed to obtain a full numerical solution, and as such it is often preferable, 
as long as we have a measure for the errors. Therefore, in order to judge whether such 
approximate solutions can reasonably be used in real applications, we have numerically 
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solved the underlying model and we have calculated errors between the approximate and 
the numerical solutions. A common characteristic of all the solutions is that the errors 
decrease with increasing axial injection velocity u c . For the PK [lj and AF [3] solutions, 
the errors are smaller than 1% for u c > 20. The KT solution yields larger errors. Each 
solution is characterized by a given profile for the density and velocity at injection, and 
the user has little freedom for adjusting the profiles for more general situations. For 
this reason, we have introduced a new class of solutions, that includes the PK and the 
KT solutions as special cases. We have shown that this new family of solutions can 
describe fairly general injection profiles, yielding errors comparable to the PK and AF 
models. Finally, we have shown that a fairly simple extension to the case of a non-null 
azimuthal velocity uq is possible. Interestingly, the effect of the azimuthal velocity is to 
distort the geometry of the plume in such a way that the density decreases more sharply 
in the vicinity of the injection boundary. In real applications, such a scenario might be 
more favorable for plasma contactors, where it is preferable to have the plume as close 
as possible to its source, in order to maximize the efficiency of the charge neutralization 
process. 
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Appendix: Derivation of Equations (11311 - (I18j) 

Substituting Eqs. (l 8 jl- (TlOjl in the continuity equation (jT]) yields: 


2 n c u c a' + a[u c n' c + n c u' c ) = 0 , 


(57) 


which reduces to 


(a 2 n c u c )' - 0 


(58) 


from which Eq. (fT3ll follows. For what concerns the momentum equation in the r 
direction ([ 2 ]) one gets: 


7 {n c n t y i n' t 
n c n 2 t 

that can be separated as 


+ pau c u 2 (a'u' c + u c a") = 0 , 


(59) 



( 60 ) 
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from which Eqs. (TIBI) and (TT5|) follow. 

Substituting Eqs. (j8|)- (TT0|) in the momentum equation in the z direction (J3]) one obtains: 


l{n c n t ) 


7-2 


n t n - pn c n t 

a 


+ u c u 2 u' c = 0 , 


(61) 


which is unfortunately not separable, and this is the source of the over determination 
of the system of Eqs. (TT3l) - (fT8l) . In order to make Eq. dfiTT) separable, one has to make 
the assumption 


'ynf 2 n c + u c u' c = 0 


(62) 


that can be integrated to obtain Equation (TT4l) . Finally, using this assumption and 
taking the p derivative of Eq. dBTj) . the result separates as 


= _ 2n Wt = D 

n c a! 2 n t u[ — (7 — 1 )n' t u t 


( 63 ) 


that ultimately yields, by integration, Eqs. (TIB]) and (TT71) . 
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Figure 1 . PK solution. Profiles for the density n and the radial velocity u r (rescaled 
by 100) at the injection boundary z = 0. In the region r > 50 where the density 
becomes imaginary, it is artificially set equal to zero. 
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u =2 u =10 



Figure 2. PK solution. Two-dimensional profiles of the density in (z,r), for 
u c = 2,10, 20,100, in logarithmic scale. The density is artificially set equal to zero 
in the regions where it becomes imaginary. 
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Figure 3. PK solution. Errors e r (black) and e z (red) as function of u c . 



Figure 4. AF solution. Profile of the density n at the injection boundary z = 0, 
for u c = 2,10,20,100 (in black, red, blue, and magenta, respectively). In the region 
r > 50 where the density becomes imaginary, it is artificially set equal to zero. 
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Figure 5. AF solution. Profile of the axial velocity u z at the injection boundary 
z = 0, for u c = 2,10, 20,100 (in black, red, blue, and magenta, respectively). Vertical 
axis in logarithmic scale. 



Figure 6. AF solution. Profile of the radial velocity u r at the injection boundary 
z = 0, for u c = 2,10, 20,100 (in black, red, blue, and magenta, respectively). 
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Figure 7. AF solution. Two-dimensional profiles of the density in (z,r), for 
u c = 2,10, 20,100, in logarithmic scale. The density is artificially set equal to zero 
in the regions where it becomes imaginary. 
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Figure 8. AF solution. Errors e r (black) and s z (red) as function of u c . 



Figure 9. KT solution. Profiles for the density n (rescaled by 2) and the radial 
velocity u r at the injection boundary z = 0. 
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Figure 10. KT solution. Profile for the axial velocity u z at the injection boundary 
z = 0, for u c = 100,10,1 (in black, red, and blue, respectively). The vertical axis is in 
logarithmic scale. 
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Figure 11. KT solution. Two-dimensional profiles of the density in ( z,r ), for 
u c = 2,10, 20,100, in logarithmic scale. 
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Figure 12. PK solution. Errors e r (black) and e z (red) as function of u c . 



Figure 13. PK and AF solutions. Angle (in degrees), with respect to the axis r = 0, 
that determines the (approximately straight) isocontour n = 0, as a function of the 
injection velocity u c . The expansion is bounded within this cone angle. 
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Figure 14. Generalized solution. Density profile at the injection boundary z = 0 for 
D = —7, -3, 1, 5 (in green, black, red, and blue, respectively). 



Figure 15. Generalized solution. Radial velocity profile at the injection boundary 
z = 0 for D = —7, -3, 1, 5 (in green, black, red, and blue, respectively). 
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Figure 16. Generalized solution. Axial velocity profile at the injection boundary 
z = 0 for D = —7, -3, 1, 5 (in green, black, red, and blue, respectively). 
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Figure 17. Generalized solution. Two-dimensional profiles of the density in (z,r), for 
u c = 2,10, 20,100, and D = —5, in logarithmic scale. 
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Figure 18. Generalized solution. Two-dimensional profiles of the density in (z,r), for 
D = —7, -3, 1, 5, and u c = 20, in logarithmic scale. 
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Figure 19. Generalized solution. Error s z as a function of uc for D = —7, -3, 5 (in 
blue, green, red, respectively), and F = 1. 



Figure 20. Generalized solution. Error e r as a function of uc for D = —7, -3, 5 (in 
blue, green, red, respectively), and F = 1. 
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Figure 21. Ratio v/ug as a function of r at 2 = 0, for E = 250 (black), 500(red), and 
1000 (blue). 




Approximate semi-analytical solutions for the steady-state expansion of a contactor plasmaAl 


E=0 E=500 



.5 

.5 

.5 

.5 

.5 

.5 


Figure 22. Generalized solution. Two-dimensional profiles of the density in (z, r ), for 
ug A 0, D = —10, F = 1, and E = 0, 500, 1000, 2000, in logarithmic scale. 
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Figure 23. Generalized solution. Two-dimensional profiles of the density in (z, r ), for 
ug A 0, D = 3, F = 1, and E = 0, 500, 1000, 2000, in logarithmic scale. 



















